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ABSTRACT 



Accepted for publication in the Astrophysical Journal 
We estimate the intrinsic neutral gas density in Damped Lyman a systems 
(n^Hi^^^) in the redshift range 2.2 < z < 5 from the DLA SDSS DR_3 sample of 
optically selected quasars. We take into account self-consistently the obscuration 
on background quasars due to the dust present in Damped Lyman a systems. 
We model the column density and redshift distribution of these systems by using 
both a non-parametric and a parametric approach. Under conservative assump- 
tions on the dust content of Damped Lyman a systems, we show that selection 
effects lead to underestimating the intrinsic neutral gas density by at least 15% 
with respect to the observed neutral gas density. Over the redshift range [2.2; 5.5] 
we find fi^^"^^ = 0.97^Q;oglQ;i5 ■ 10~^, where the first set of error bars gives the 
1(7 random errors and the second set gives the modeling uncertainty dependent 
on the fraction of metals in dust - from 0% to 50%. This value compares with 
^^Hi"^^ = 0.82^gQ5 (la error bars), which is obtained when no correction for 
dust is introduced. In the model with half of the metals mass in dust we cannot 
constraint at a confidence level higher than 90%. In this case there is 

indeed a probability of about 10% that the intrinsic column density distribution 
of DLA systems is a power law /(A///) oc In contrast, with 25% of the 

metals in dust - the most realistic estimate - a power law is ruled out at 99.5% 
of confidence level. 



Subject headings: dust, extinction - galaxies: high- redshift - intergalactic medium 
- galaxies: ISM 
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1. Introduction 

Damped Lyman a systems (hereafter DLA systems) are quasar absorption systems 
with a column density above 2 ■ 10^°cm^^ and represent the high end of the distribution of 
absorption systems starting from the Lyman a forest at Nhi ^ lO^'^cm'^. DLA systems 
represent the most significant reservoir of neutral hydrogen in the universe available for star 
formation. These systems are considered to be either (cold) massive rotating disks, the 
progenitors of todays disk galaxies (Prochaska & Wolfe 1997; Wolfe et al. 2005) or compact 
protogalactic clumps (Haehnelt et al. 1998; Nagamine et al. 2004). 

In the era of precision cosmology, an accurate measure of the total mass density of neu- 
tral gas as a function of the redshift represents an important constraint for galaxy formation 
models. Ground observations are able to identify DLA absorption features in the spectra 
of quasars from Zabs ^ 1-8, where the absorption lines enter the atmospheric window, up 
to Zabs ~ 5.5. With an all sky survey like the Sloan Digital Sky Survey, spectra of several 
thousands of quasars with enough resolution for DLA detection have been acquired and the 
observed density of neutral gas in DLA systems is now measured with errors below 10% 
(Prochaska et al. 2005). 

This measurement must be interpreted with some caution, as the presence of DLA 
systems along a line of sight leads to a potential obscuration due to the dust that they 
host: the observed gas density is a biased estimator of the intrinsic density unless the dust 
effects are accurately quantified. Several papers, starting from Ostriker & Heisler (1984) 
have attempted to model the influence of dust along the line of sight, often with conflicting 
results. 

A detailed analysis framework for the obscuration of quasars has been developed by Fall 
& Pei (1993) (see also FaU & Pei 1989; FaU et al. 1989; Pei et al. 1991; Pei & FaU 1995) and 
applied to the quasars sample of Lanzetta et al. (1991). Their study highhghted a potentially 
severe effect of the dust bias that did not allow to put an upper limit to the intrinsic density 
of DLA systems. In fact, absorbers with high column densities and/or with high dust-to-gas 
ratio represent essentially "bricks" along the line of sight to a quasar and are very likely 
to be missed in optically selected surveys. An additional evidence for the dust obscuration 
came in the form of a detected preferential reddening in the spectra of quasars with DLA 
absorption with respect to a control sample without detection of these systems (Pei et al. 
1991). 

More recent investigations revised these earlier results on the dust content in DLA 
systems and on their reddening of background objects (Murphy & Liske 2004; Ellison et al. 
2005), finding in particular no robust evidence for the reddening of quasars at 2; ~ 3 with 
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DLA features in their spectra: at 3a Murphy & Liske (2004) find E{B-V) < 0.02mag, while 
EUison et al. (2005) have E{B — V) < O.OAmag . At the same time radio selected quasars 
surveys (Ellison et al. 2001), with complete optical follow-up detection have provided the 
first bias free constraints on the intrinsic distribution of DLA systems. 

Taking advantage of these recent measurements for the number density of DLA systems, 
we have previously characterized (Trenti & Stiavelli 2006) the dust absorption along random 
lines of sight by means of a Monte Carlo code, finding that, on average, the deviations from 
unit transmission are effectively modest ((exp(— r)) > 0.9 at an emitted wavelength Ae = 
0.14/7,771 over all the redshift range) and of limited impact on most observations. However, 
this result does not exclude the presence of a small fraction of lines of sight (of the order 
of a few percent) through the most massive and/or the most metal-rich DLA systems and 
characterized by a large optical depth. Indeed, Wild & Hewett (2005) and Wild et al. 
(2006) find a significant evidence of reddening in DLA systems with Call absorption fines 
at moderate redshift {zabs ~ 1; {E{B — V)) > 0.1). Similarly York et al. (2006) measure 
E(B — V) up to 0.085 for Mgll selected DLA systems at Zabs ~ 1-0. As the determination of 
the gas density in DLA systems is dominated by these most massive absorbers, the potential 
bias in this measure, correctly stressed by Fall & Pei (1993), must not be dismissed by the 
recent evidence of a very modest average deviation from unity transmission. 

In this paper we take advantage of the large sample of DLA systems identified in Sloan 
quasars (Prochaska et al. 2005) and we investigate the relation between the observed and 
intrinsic density of neutral gas in these systems. The Sloan sample that we consider has 
three main advantages over the sample used by Fall & Pei (1993). (1) It is larger by a factor 
30. (2) The quasars have been selected within a luminosity limit in the I band, significantly 
less sensitive to dust obscuration than the B band used for the older sample. (3) The color 
selection algorithm has a good sensitivity to red quasars (Richards et al. 2002). In fact 
even if a quasar with a DLA absorption system is not dropped below the I-band fiux limit, 
its color can be changed so that it ends up lying outside the color selection box used to 
identify quasar candidates for follow-up spectroscopy. Thanks to the precision of the Sloan 
photometry, that allows a clear separation of the stellar locus in the color space, and to the 
use of an extended color selection box, the loss of completeness due to this latter effect is 
only marginal (see Richards et al. 2002 for a detailed discussion of the Sloan color selection 
algorithm and for completeness tests) and is therefore not considered in this work. 

This paper is organized as follows. In Sec. 2 we characterize the obscuration bias in a 
magnitude limited survey, in Sec. 3 we describe the dataset that we are using. In Sec. 4 we 
present our analysis for the parametric estimation of the intrinsic comoving density of neutral 
gas, whose uncertainties are quantified in Sees. 5-6 by means of Monte Carlo simulations of 
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synthetic observations. In Sec. 7 we discuss the accuracy of a non-parametric estimator for 
the neutral gas density. We summarize our findings in Sec. 8. 



2. Dust effects in magnitude limited surveys 

In a magnitude limited survey, obscuration along the line of sight leads to the potential 
loss of some hues of sight. Following Fall & Pei (1993), to quantify the effect let us consider 
an intrinsic quasar luminosity function ^i{L) in the interval [L^i„, +00], where is the 
luminosity above which the sample is complete. The intrinsic number of objects above the 
completeness threshold is: 

r+00 

Ni= / ^i{L)dL. (1) 



The presence of dust introduces an optical depth r along the line of sight, with a corre- 
sponding transmission coefficient exp{— r}. Under this condition the number of objects 
with observed luminosity above L^i„ is: 

r+00 

N,= ^,{L)dL. (2) 

The dust obscuration leads to a fraction of missing objects given by No/Ni. 

No lLl.ML)dL 



00 



(3) 



In general the precise value of No/Ni for a given r depends on the sensitivity limit of the 
survey Lmin and on the form of the luminosity function ^{L). Under the assumption that 
the luminosity function is a power law $(L) oc L^^^^ (in the luminosity range from L^j^ to 
+00) the ratio No/Ni can be easily computed and is independent of Lmin'- 

Ni~ l:' ~ ^ ' 



mm 



This was already noted by Fall & Pei (1993). 

The luminosity function for quasars is modeled in terms of a gamma function and/or 
of a double power law (e.g., see Pei 1995). However the SDSS spectroscopic quasar sample 
is shallower, for Zqso ^2.2, than the knee of the distribution so that the luminosity function 
for the quasars that we are considering in this paper can be effectively treated up to L^in 
as a single power law (see Richards et al. 2006), simplifying significantly our analysis by use 
of Eq. (4). 



- 5 - 



In the following analysis we assume a redshift dependent luminosity function (Richards 
et al. 2006): 

Mh,z)^A{z)Lf^'^-\ (5) 

where L/ is the luminosity in I band^ (where the sensitivity limit for detection of quasars 
in SDSS is given) and f3{z) is a slowly evolving function of the redshift with /3 varying from 
2.2 to 1.1 in the redshift range [2; 5] (Richards et al. 2006): 

I3{z) = 2.1 - 0.275(z - 2A5)e(z - 2.45), (6) 

where 9{x) is the step function defined as 9{x) = 1 if x > and 9{x) — otherwise. The 
SDSS DLA DR_3 sample of quasars considered in this paper has an average redshift of 2.97 
with a standard deviation of 0.65. The average slope can be approximated as: 

{(3) = 1.95 ^ 2. (7) 



In case of absorption due to a dusty DLA system at redshift Za with column density 
Nffi, the optical depth at an observed wavelength Xg can be written as: 

t{Xo, Za, NhI, k) = kNHli{Xol{l + Za)), (8) 

where k is the dust-to-gas ratio and ^ the relative extinction curve normalized to the ab- 
sorption cross section in B band cr(As) (e.g., see Pei 1992): 

e(A) = a{X)/a{XB). (9) 

If Nhi is expressed in units of lO^^cm"^, the absorption cross section in B band is kNni with 
k — 0.8 for galactic dust. The value of k for DLA systems depends on their metallicity Z and 
on the fraction of metals in dust, i.e. on the dust-to-metals ratio. DLA systems are generally 
characterized by a low metallicity and by a moderate evolution of their properties with the 
redshift (Wolfe ct al. 2005). We approximate the observed average (Hl-column density 
weighted) rcdshift-mctallicity relation (based on a number of observations, e.g. Prochaska et 
al. 2003, Dessauges-Zavadsky et al. 2004, 2006, Akerman et al. 2005) with a linear function 
for \og[Zo{z)). This provides a good agreement with the data (e.g. see Fig. 13 in the 
compilation by Kulkarni et al. 2005) in the redshift range 2 < 2; < 5: 

Z<,(2)/Zq = 0.2- 10-°-2^ (10) 



^in this paper we assume an average wavelength for the band of (A/) = O.S/xm. 
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The average metallicity in Eq. (10) translates, for a Milky Way dust-to-metals ratio (i.e. 
50% of the metals in dust grains), into a "Milky Way" average dust-to-gas ratio: 

W(^) =0.16-10-°-'^ (11) 

To account for different fractions of metals in the form of dust grains we introduce a correction 
factor for the intrinsic dust-to-gas ratio ki{z): 

ki{z) = ai^kuwiz). (12) 

DLA systems are considered to have a smaller fraction of metals in dust than our galaxy 
(with roughly one quarter of the total metals content in dust) , so realistic values for are 
expected to be around = 0.5 (Pettini et al. 1997; Vladilo 2002; but see Pei et al. 1999 
where the Milky Way dust-to-metals ratio is used). In the following sections we present 
our analysis using a range of dust-to-metals ratio, highlighting the dependence of the dust 
bias on this quantity. We assume a reference value = 0.5 (25% of metals in dust, the 
value derived by Pettini et al. 1997), but we include a grid of models with G [0; 1] with 
Aa^ — 0.125, discussing in depth also the case = 1 (50% of metals in dust). 

Our analysis is carried out by assuming that all DLA at a given redshift have a fixed 
metallicity, as well as a fixed dust-to-metals ratio. A detailed modeling of the metallicity and 
dust-to-metals ratio distributions for the purpose of determining Vf^l"^^ is extremely chal- 
lenging, as only a small subsample of DLA system have measured metallicities. Fortunately, 
the value of fi^^'^'' is not expected to significantly depend on the scatter in the dust-to-gas 
ratio ki of the absorbers (see Fall & Pei 1993, Appendix A). To ensure that this is the case, 
in Sec. 6 we validate our scatter-less approach by analyzing synthetic observations with a 
variety of dust-to-gas distributions. 



3. The data: The DLA SDSS DR_3 and the CORALS surveys 

The data used in this paper are primarily taken from the DLA survey by Prochaska et 
al. (2005). We consider all the quasars from their Table 1 and all the DLA systems reported 
in their Table 3. The sample consists of 525 DLA system found in 4568 spectra of quasars 
with a minimum signal to noise ratio of 4. 

In addition we consider the DLA systems detections in the radio selected CORALS 
survey (EUison et al. 2001), intrinsically free from dust bias. This survey consists of 19 DLA 
systems detected in 66 spectra of quasars (see their Table 3) . The statistical sample for the 
survey (see EUison et al. 2001), that we consider for the analysis, is restricted to 17 DLA in 
the redshift range [1.8; 3.5]. 
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4. Parametric Estimation of Q^^^^ from SDSS DR3 

Following the standard practice (e.g. Lanzetta et al. 1991) we define the number of 
observed DLA systems in the intervals [Nhi, Nhi + dNni] and [X, X + dX]: 

foiNHi,X)dNHidX, (13) 

where fo is the observed frequency distribution and dX the absorption distance: 

dX^^^il + zfdz. (14) 

In this paper we adopt the standard KCDM cosmology with Hq — 70km/s/Mpc, Q\ — 0.7 
and Qm — 0.3, so that: 

dX = , ^ ^ =dz. (15) 
V0.7 + 0.3(1 + ^)3 

The data from SDSS DLA DR_3 span over a total integrated absorption path-length AX — 
7333.1. 

Considering the obscuration bias discussed in the previous section for the Sloan survey 
(see also Fall & Pei 1993) and assuming that the distribution of absorbers along the line 
of sight is not correlated^, wc can write the relation between the intrinsic and the observed 
frequency distribution of DLA systems as: 

f„{NHi, z) = MNhi, X(z))e-<^>"'^'='^^(^)^«^^(^^/(^+^)), (16) 

where we assume the average slope {f3) = 2 for the quasar luminosity function. Following 
Prochaska et al. (2005) we parameterize / using a gamma function: 

,(^,„X)^,(^)-exp(-^), (IT) 

and a power law: 

MNHi,X)=7],iNHir'. (18) 

The use of the gamma function for the intrinsic distribution /j has the advantage that, given 
Eq. (16), also the observed distribution remains a gamma function. A intrinsic power law 
is instead mapped into a observed gamma function by the effect of dust. 



^We are essentially neglecting clustering along the line of sight. 
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If we assume that the intrinsic distribution of DLA systems follows a Poisson statistics, 
the observed distribution will also follow a Poisson statistics. The likelihood to be maximized 
is therefore (see also Eq. 27 in Fall & Pei 1993): 

dz dNni ■ UNhis, ^s)-r\Il foiNmn, Zn), (19) 

where the sum extends over all the sample of quasars that have been searched for DLA 

systems in the redshift interval [zmf^', Zgupg], where the signal to noise ratio is above 4, and 
the product is over all the DLA systems detected in those quasars. In Eq. (19) the observed 
distribution function is expressed in terms of the intrinsic distribution /j through Eq. (16). 
Ndla = 2 • 10^^ cm"^ is the assumed lower limit for the column density of a DLA system. 

Before presenting the results from the maximum likelihood analysis we investigate the 
properties of with the aim to discuss upper limits on the dust-to-gas ratio and on Vt^^['^\ 



4.1. Basic considerations from the properties of 

Fig. 1 shows the observed column density distribution of DLA systems, averaged over 
redshift. The observed data are well fitted by a gamma function (Prochaska et al. 2005, see 
also Eq. 17 and Tab. 1) with a knee at (iV-y)o — 3 • lO^^cm"^ and a slope (mi)o — —1.8. 

The origin of the observed knee at {N^)o can in principle be either due to an intrinsic 

decrease of the distribution at high column densities, i.e. due to iV^, or due to the effect of 
obscuration bias, that introduces in the data an exponential-like decrease. We can obtain a 
limit on the dust-to-gas (expressed in terms of ««) ratio assuming that the observed knee is 
entirely due to dust bias. From Eq. (16) we can write: 

{(5)a,kMw{z)i{\i/{l + z)) < l/{N^)o, (20) 

so that, considering that the average redshift of DLA absorbers in our sample is {z) = 3.1 
we estimate a limit on a,^: 

^1 1 ^1 1 ^ 

^ WyVo ■ mMw{{z)mi/{l + {z))) - 3 ■ 2 • 0.038 • 2.47 ~ ^ 

As this value implies about 90% of metals in dust, which appears extremely improbable for 
DLA systems, the presence of an observed knee in the column density distribution of the 
gas is likely an intrinsic feature and not induced by the dust bias alone. This implies that 
the expected difference between the intrinsic and the observed density of neutral gas in DLA 
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systems is limited for the most realistic fraction of metals in dust (25% i.e. = 0.5). Only 
in the case of a larger fraction of metals in dust (i.e. ~ 1) there can be the possibility of 
a significant dust bias, as in that case an acceptable model of the intrinsic distribution can 
still be obtained in terms of a power law. In the next Section we quantify more precisely the 
dust bias by studying the likelihood L for /j at increasing ^^^^^ for different values of 



4.2. Whole sample analysis 

By aggregating all the data in SDSS DLA DR_3, the results of our maximum likelihood 
analysis are reported in Tables 1-2 and in Fig. 2. The results for the fit with an intrinsic 
gamma function are reported in terms of the comoving density of gas in DLA systems 
^m^^\ that is the first moment of /«: 

^^BLA)^^ = ^f!^ dNHiNHiMNni, z)dX, (22) 

^Pc JNdla 

where ruH is the mass of the hydrogen atom, fi — 1.3 is a correction factor for the composition 
of the gas, c the speed of hght, Pc the critical density and Ndla — 2 • 10^° cm~^. 

We obtain the following results: 

• = 0: If we neglect the effect of dust absorption so that fo = /i, we re-derive 
the parameters for the gamma distribution reported in Prochaska et al. (2005) (see 
our Tab. 1). All our likehhood values for the Sloan data are presented in log units 
normalized to the likehhood value for the best fitting dust-free model. For — Q 
the likelihood is sharply peaked around the best value for fi^^^-*. However, when we 
introduce the effect of the dust («« > 0) there is an increasingly strong degeneration 
toward high values for O^^'^'* (see Fig. 2). 

• a^ — 0.5: For our standard scenario the likelihood has a maximum at Q^^"^^ = 0.97 • 
10~^, to be compared with the value 

(^i??^^)^ = 0.82 • 10-^ given by the maximum 
likelihood analysis with a^, — Q (no dust). In our modeling of the column density 
distribution the metallicity is introduced has an external input, based on independent 
observations. An additional support to motivate the need of taking into account the 
obscuration bias comes from the fact that the likelihood for = 0.5 is higher than 
for = 0. 



• ckk = 1: The bias induced by the dust increases as larger are considered. As 
expected qualitatively by the simple considerations on the shape of fo presented in 
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the previous section, the maximum likelihood analysis for our large dust-to-metal ratio 
scenario = 1) is unable to set a strong upper limit to fi^^"^^ as there is a wide 
wing of high likelihood values toward large neutral gas densities. At = 1 a power 
law fi with slope m2 = —1.95 has a likelihood ratio over the best gamma solution 
of 7^ = 2\og{Cr/ jC-n^) = 2.1. A likelihood ratio test with this value implies that a 
power law solution can be ruled out only at about 85% of confidence level if we assume 
that TZ is distributed as a with one degree of freedom. However this is in general 
valid only when the model is linear in the parameters and the errors follow a gaussian 
distribution (Lupton 1993). As our model is strongly non linear, in the next section 
we estimate the confidence level by means of Monte Carlo simulations. 

• > I- If we consider fits for at q;«; > (see Table 2) we obtain an increasingly bet- 
ter modeling in terms of power laws up to ~ 2. As expected from the preliminary 
examination of /<,, eventually the likelihood decreases at higher dust-to- gas ratio start- 
ing from > 2. We recall that > 2 is formally unphysical, as it implies more than 
100% of the metals in dust, however, due to the fact that we have adopted a simple 
fit to the observed metallicity, a high value also means that the intrinsic metallicity 
content of DLA system is underestimated in our model. As it appears unlikely that 
the metallicity measurements do significantly underestimate the dust content of DLA 
systems, we expect that realistic values of < 1- 

• CORALS: For comparison we have repeated the maximum likelihood analysis on the 
radio selected CORALS quasars (Ellison et al. 2001). The results are reported in 
Table 3 and in Fig. 3. For these data the measure is guaranteed to be free from dust 
bias, but the sample size is insufficient to accurately constraint the parameters of 
(the likelihood curve in Fig. 3 is relatively flat around the maximum). The maximum 
is at n^^j^^^ = 1.2 • 10~^ with extremely large la uncertainties (of the order of 50%) 
due to the flatness of the likelihood function. In addition, the maximum likelihood 
for the gamma function model is only negligibly better than the one for a power law 
(7^R^ 0.1). 



5. Errors estimate 

To estimate the error on the determination of Q^^j^"^^ we have performed a bootstrapping 
analysis of the data. We have randomly extracted samples of lines of sight with uniform 
probability from the SDSS DLA survey. The samples have the same number of lines as the 
original data. For each of these simulated samples we have performed a maximum likelihood 
analysis to determine In the case oi a,^ — 1, the results obtained with 400 random 
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realizations are presented in Fig. 4 and show that Q^^"^-* is indeed distributed around the 
best value derived from the original data (1.25- 10^'^). which coincides with the median of the 
distribution. The average value of the distribution is 1.27- 10^^. At 68% of confidence level 
the measurements lie in the interval [1.10 ■ 10~^; 1.41 ■ 10~^]. We assume this interval as our 
fiducial error on r^if/^^. At 90% of confidence level (^Jf/^^ < 1.48 -10^^ The bootstrapping 
procedure has been repeated for all the scenarios with different a^, resorting to at least 100 
simulations for each value of the dust-to-metals ratio considered. The results are reported as 
standard la error in Table. 1. The case = 0.5 is shown in Fig. 5 (using 200 simulations 
in this case). 

As a additional test to investigate the uncertainties in our analysis, we have simulated 
a number of synthetic samples of data using the Monte Carlo code developed in Trenti & 
Stiavelli (2006). These samples have then been processed using the same method apphed 
to the real observed data, to test its efficiency in retrieving the input parameters and to 
quantify the typical error. 

Each random realization of a simulated observation is characterized by the same ob- 
served path- length of the SDSS sample {AX — 7333), but consists of lines of sight with 
DLA distribution (and detection) in the redshift range [2.2; 4]^ The optical depth along each 
intrinsic line of sight is computed from the realized DLA distribution and then a test on the 
obscuration probability (Eq. 4) is performed to accept or reject the intrinsic line of sight as 
a observed one. The only source of uncertainty in these synthetic observations comes from 
the discrete sampling of fi and from the obscuration probability test. The model does not 
include observational errors. 

We have simulated two models: (i) an intrinsic gamma function and (ii) an intrinsic 
power law, with the parameters given by the maximum likelihood analysis from the SDSS 
data (see Table 1 for the gamma function and Table 2 for the power law). 

Fig. 4 shows, for = 1. the distribution of Q^^"^^ as measured from the simulated 
observations (80 different realizations) of a intrinsic gamma function. It is reassuring that 
the distribution basically agrees with the one obtained using the bootstrap method. The 
mean value for il^^^-* is 1.29 • 10~^ with a standard deviation of 0.23 • 10^^. As we used an 
input value Q^'^'^^ = 1.25 • 10~^ to generate the synthetic observations, the fitting procedure 
is able to recover the intrinsic density of neutral gas introducing only a negligible bias, which 
is much smaller than the one sigma random uncertainty. 



^This has been done for computational reason in order to speed up the evaluation of the double integral 
in Eq. 19. 
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The analysis of the simulated observations generated from an intrinsic power law distri- 
bution allows to quantify the probability that the identification of a gamma function as best 
fitting model for the Sloan data is only due to the combined effect of discrete sampling and of 
the additional degree of freedom of the gamma function over a power law. For our reference 
scenario {a^ = 0.5) we have generated 400 random realizations of observations starting from 
an intrinsic power law with unconstrained neutral gas density (772 = 1.54- 10~^, 1712 = —2.05). 
The likelihood ratio (gamma over power law) distribution shows only 2 realizations that have 
TZ greater than the one measured from the SDSS DLA data. This means that an intrinsic 
power law can be ruled out at 99.5% of confidence level. Repeating the experiment in the 
case of — 1, r]2 — 1.77 • 10~^ and m2 = —1.95 we find that about 8 % of the realizations 
have a value of TZ greater than the one measured from the SDSS DLA data. This means that 
the confidence level at which a power law distribution (with unconstrained upper cut-off) 
can be ruled out is in this case 92%. 



6. Synthetic observations with dust-to-gas ratio scatter 

In the maximum likelihood applied in this work we assume that all the absorbers have 
a fixed dust-to- metals ratio a^^ and a metallicity dependent only on the redshift. In reality a 
significant scatter around the mean values is expected for both these quantities and has in- 
deed been observed (e.g. see Kulkarni et al. 2005). Fall & Pei (1993) have already noted that 
the scatter in the dust-to-gas ratio is not introducing significant bias under the assumption 
that there is no correlation between fcj and the column density Nhi- 

To confirm the results of Fall & Pei (1993) and to extend the analysis of the bias to 
the case where a correlation between ki and Nhj is present we have simulated a number of 
synthetic samples of data adopting an approach similar to the one described in Sec. 5. In 
this case we compute the optical depth of each absorber by drawing dust-to-gas ratio ki from 
a given distribution, rather than adopting a fixed value depending only on z. 

When there is no correlation between fcj and Nhj we use (based on Eq. 12): 

ki{z) ^ ^ ■ a^- kMw{z): (23) 

where ^ is a random variable with a given probability distribution p{C,) such that the expected 
value of ^ under p is Ep[^] = 1. We have considered the following form for p(^): 

• A uniform distribution around the average value: p(^) = 1 if ^ G [0.5; 1.5] andp(^) = 
otherwise. 
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• A lognormal distribution with variance parameter cr*, as studied by Fall & Pei (1993): 
piO = exp {- log2 (eexp {a2/2})/(2a2)}/y2^^. Wc adopt a, = 0.5, 1, 1.5 inspired 
by the dispersion in the dust-to-gas ratio observed for galaxies in the local universe 
(Pei 1992). 

In addition wc have generated synthetic observations with a correlation between fcj and Zi 
adopting the following prescription: 

JNdla ^^^^^ ■ ^^HI JiV^Hl) 

where 7 = —0.8, —0.4, 0.4, —0.8 (i.e we explore both a increasing and decreasing trend of 
ki vs. Nhi) and ^ is a random variable with the same properties bescribed above. The 
normalization adopted is such that the HI column density averaged value of ki is o;^ ■ kj^^/yzj . 

For each combination of the parameters considered (see Table 4), we have generated 
60 synthetic SDSS-like samples of data (as in Sec. 5) adopting a gamma function model for 
fiNni) {N^ = 4.15 ■ lO^^cm"^; mi = —1.8) and = 0.5. These synthetic samples have all 
Qjf/^) = 0.97 • 10-^ By applying our maximum likelihood we obtain the following results 
(summarized in Table 4): 



• Uncorrelated scatter in ki. The maximum likelihood analysis is able to correctly recover 
the value of fi^^'^'' used to generate the observations (see Fig. 6 and Table 4). There 
is only a small bias (at the percent level) toward obtaining output values marginally 
larger than the input, even when the dust obscuration is such that the observed column 
density averaged dust-to-gas ratio is smaller by 15% with respect to the intrinsic value, 
as happens in our lognormal model with a* = 1. This demonstrates that our analysis 
in terms of the average dust-to-gas ratio is robust with respect to uncorrelated scatter 
in ki- Interestingly we also note that the difference between the observed and the 
intrinsic dust-to-gas ratio for the lognormal model with cr* = 1 is in agreement with 
the difference between the metallicity in optically and radio selected samples of DLA 
systems (see Akerman et al. 2005)^. 

• Correlated scatter in ki. In this scenario the analysis in terms of an average dust-to-gas 
ratio systematically mis-estimates the dust content at both ends of the column density 
distributions of the absorbers. In particular, if ki increases with Ni, then absorbers 



^Note however that the difference in the metalUcities measured by Akerman et al. (2005) is not particularly 
robust from a statistical point of view, as it is within the la error bar. 
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with the largest column densities will be more dust rich than assumed by our analysis. 
This means that these absorbers will more likely be dropped by the sample than is 
assumed in the modeling, so that is underestimated. When ki decreases with 

Ni, the opposite happens: the densest absorbers are less dusty than assumed and 
fl^^j^^^ is overestimated. Due to the skewness of the lognormal distribution the bias 
introduced by the correlation is expected to be larger when ki is a decreasing function of 
Nhi- This is indeed what we observe in the analysis of our synthetic observations (see 
Fig. 7 and Table 4). The bias introduced is, however, relatively modest (at most 15% in 
^i?/^^^)' ^^^^ assuming a rather strong dependence of ki on Nhi (i.e. ki oc N^f^). The 
bias could be avoided by taking into account the precise form of the correlation into 
the modeling of the the dust-to-gas ratio. Unfortunately the ki vs. Nhi correlation 
appears hard to quantify with the present data. Only a small subset of the known 
DLA systems have a measured metallicity and highly dusty and dense absorbers are 
preferentially missed from the sample. In addition there arc some evidences that the 
fraction of metals depleted in dust depends on the column density of the absorbers 
(Welty et al. 1997). Therefore plots of the measured gas phase metallicity versus the 
hydrogen column density, such as those presented in Boisse et al. (1998) and Savaglio 
(2001) cannot be easily interpreted to extract the ki{NHi) relation. 



7. Non-parametric Estimation of ojf/^^ from SDSS DR3 

The comoving density of DLA systems can be also estimated using a non-parametric 
approach by taking the discrete hmit of Eq. (22) combined with Eq. (16): 

{AJq^- ^^^^ 

The sum is done over all the DLA systems in the sample, while the intrinsic total path-length 
{AX)i is given by summing over all the individual path-lengths for the quasars in the survey, 
with a weight depending on the optical depth r„ along the line of sight to the n — th quasar 
and on the redshift-dependent slope of the quasar luminosity function I3{z): 

(AX), = Y.^dXo)ne^-^- (26) 

n 

The optical depth r„ is given by: 

Tn = ^ kMw{Zp)NHIpC{Xl/{l + Zp)), (27) 

p 

where the index p is running over all the DLA systems along the n — th line of sight. 



- 15 - 



The random uncertainty on the measure has been estimated by using a 'Jack-knife' 
analysis (Lupton 1993): we derive fi^^'^'' for N = 22 subsamples of quasars by ignoring 
Ntot/N quasars each time. The uncertainty on O^^'^'' is then given by: 



^' = ^^E(^^- Wf, (28) 



N 



where stands for Ol?/"^'*„. 



HI 

We obtain the following results: 

» = 1: Eq. 25, apphed to the whole sample with = 1, gives Q^^^^^ — 1.04 • 10""^ 
with an error of 8 • 10~^ (see Table 5). The obscuration bias based on this estimate is of 
about 30 %. This represents a lower limit on the bias, as from Fig. 8 it is evident that 
with = 1 the dust bias is so severe that the discrete evaluation of Q^^^"* has not 
converged. On the top of the random errors, estimated with the 'Jack-knife' method, 
this non parametric estimation suffers from systematic effects of the order of at least 
10 - 20%. 

» «K = 0.5: The values of Q^^"^^ obtained in our reference model with = 0.5 for 
different redshift intervals are reported in Tab 5). The systematic corrections due to the 
dust bias arc relatively modest (of the order of 10%) but still represent a contribution 
to the total uncertainty that cannot be neglected. Like in the case discussed above a 
systematic effect is also present. Here the magnitude of the effect is smaller, of the 
order of 5%, as can be evaluated from the difference in between this method 

and the maximum likelihood analysis of the previous section. 

• CORALS: We also apply the same analysis on the DLA detections in the CORALS 
survey (Ellison et al. 2001), obtaining fi^^^-* = (1.2 ± 0.5) ■ 10^'^ (la error bars; see 
Table 5). Note that our analysis leads to a different value for fi^^"^^ because in Ellison 
et al. (2001) the analysis was performed with a different cosmology [VLm = 1, J^a = 0). 
The main limit of the survey is its modest path-length extension (AA = 200.8, to be 
compared with the Sloan pathlength AA = 7333.1), which is translated into a large 
random error on fi^j. In addition potential systematics errors due to incompleteness 
may bias the CORALS measure much like dust introduces a bias in optically selected 
surveys. In fact the likelihood modeling of the data is inconclusive, allowing a variety 
of different forms for /j (see Fig. 3). 
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8. Conclusions 

In this paper wc have improved the analysis of the column density distribution of DLA 
systems in the SDSS DR_3 DLA survey(Prochaska et al. 2005) to take into account the bias 
due to dust obscuration along the line of sight. A first modeling of the bias was constructed 
by Ostriker & Heisler (1984) and improved by Fall & Pei (1993). These earlier estimates 
expected a severe effect of obscuration, with the observed gas density of DLA systems being 
up to several times smaller than the intrinsic one. The best estimate for Q^^"^^ given by Fall 
& Pei (1993) is 4.9- 10~^ (obtained using a cosmology with flM = 1 and Hq — 70km/s/Mpc) 
and their upper limit Q^^j^^^ < 3 ■ 10~^, obtained from a standard big bang nucleosynthesis 
abundance model. More recent works (Murphy & Liske 2004; Ellison et al. 2005) tend to 
dismiss the issue of obscuration bias based on the absence of systematic reddening in the 
spectra of quasars with DLA absorption features. 

Here we show that the effect of obscuration, while not being as severe as predicted 
by Fall & Pei (1993), does indeed play an important effect on the precise measurement of 
l^if/-"). In the era of precision cosmology, where the observed density n^^^-* is constrained 
with errors below 10 %, the systematic effects are not to be underestimated. With the 
typical amount of metals present in DLA systems the observed density Vf^j'^^ derived from 
shallow magnitude limited surveys of quasars underestimates the intrinsic density Q^^"^^ 
by about 15 % assuming dust-poor DLA systems (i.e. systems with a fraction of metals in 
dust of 25%). Our best estimation for z G [2.2; 5.5] gives an intrinsic neutral gas density 
^Hi^^"* ~ 0-97^006 ■ 10~^ (lo" error bars) to be compared with the observed gas density 
i^Jf/^^ = 0.817l[|:[j^^ ■ 10-3 derived by Prochaska et al. (2005). If we leave the dust to metal 
ratio parameter free to vary over the relevant range (from to 1), we find ^2^''""^^ = 
0.97]!]Q;Qg^Q;^5 • 10"^, where the first set of error bars gives the la random errors and the 
second set gives the modeling uncertainty dependent on the fraction of metals in dust. The 
obscuration bias therefore represents the main source of uncertainties for the determination 
of the intrinsic neutral gas content in DLA systems. 

Our analysis has been carried out assuming that all the DLA absorbers at a given 
redshift have the same dust-to-gas ratio. By means of monte carlo simulations of synthetic 
observations we shoow that this assumption does not introduce a significant bias as long as 
the dust-to-gas ratio is not correlated with the hydrogen column density. A correlation of 
the form ki oc N]jj introduces a systematic error of the order +14% in Q^''^'^^ for 7 = —0.8 
and of —3% for 7 = +0.8 (the bias is reduced to +8% and 0% for 7 = +0.4 respectively). 

Caution is also needed when taking the result from the maximum likelihood at face value. 
In our reference scenario we find that a power law function for column density distribution of 
DLA systems can be ruled out at a confidence level no greater than 99.5%. If the dust content 
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in DLA systems is higher, the dust bias becomes more significant. For 50% of the metals in 
dust grains the present data do not allow to put an upper limit to the neutral gas content 
with confidence level greater than 90%. We show in fact that there is a probability of about 
10% that the SDSS DLA DR_3 data are consistent with an intrinsic power law distribution 
in column density adopting = 1. The slope of this power law is m2 = —1.95 and an 
upper cut-off cannot be derived from the data. In optically selected surveys, absorbers with 
high column densities of neutral gas and with high metallicity are missed, letting systematic 
uncertainties go easily out of control. 

Radio selected quasar samples would represent an elegant, bias free sohition to the 
measurements of ^2^^^'' . Unfortunately, with the present extension, these surveys constraint 
^ifi^^^ with a random uncertainty of the order of 40 % (plus potential systematic effects due 
to the hmited number of observed DLA systems). To reduce the errors within the 10 % level 
without the need to assume a modeling for the dust bias it is therefore necessary to increase 
the path- length of radio surveys by a factor 10 at least. The UCSD radio survey (Jorgenson 
et al. 2006) has recently made promising progresses in this direction. 

It is a pleasure to thank Michael S. Fall for stimulating discussions and Jason X. 
Prochaska for useful suggestions and comments on a draft of the paper. We are grateful 
to the referee for constructive suggestions. This work was supported in part by NASA 
JWST IDS grant NAG5-12458. 
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1 10 

N^/l02icm-2 

Fig. 1. — Observed column density distribution of DLA systems in the SDSS DLA DR_3 
sample compared to observed (solid) and intrinsic (dotted) best fitting gamma functions. 
The difference between /j and fo is due to dust obscuration, especially apparent for high 
column densities. The parameters value for the fitting functions are reported in Tables. 1-2. 
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Fig. 2. 

systems VL^j^^^^ 
= 0.5 



Maximum likelihood £ for the intrinsic comoving neutral gas density in DLA 
for the SDSS DLA DR_3 data (solid line with = 1; red dotted with 
The green dashed line is associated to the observed gas density (a^ = 0). The 



likelihood curves have been obtained by maximizing C over A^.^ and mi at fixed Q 
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and 
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Fig. 3. — CORALS survey analysis. Right: likelihood curve in function of for a 

gamma function modeling of the data. The curve shows only very limited variations around 
the maximum value. This is because the dataset is too small to allow to constraint the 
gamma function parameters. Left: maximum likelihood C for the power law slope m2. The 
likelihood is in units of the maximum likelihood value for the fit with a gamma function (see 
Table 3). 
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Fig. 4. — Q^^j^"^^ distribution for = 1 from the bootstrapping analysis of the data (red 
dotted Une; 400 reahzations) and from 80 Monte Carlo realizations (solid line) of synthetic 
observations starting from the best fitting model for /j (gamma function with mi = —1.8, 

1.25 ■ 10-3). 
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Fig. 5. — n^^"^'* distribution for = 0.5 from the bootstrapping analysis of the data (red 
dotted hne; 200 reahzations) and from 60 Monte Carlo realizations (solid line) of synthetic 
observations starting from the best fitting model for /j (gamma function with mi = —1.79, 

0.97-10-3). 
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Fig. 6 - 

of synthetic observations starting from the best fitting model for (gamma function with 
mi = -1.79, = 4.15 ■ W^^cm-^ and fi^f/^^ = 0.97 ■ IQ-^). The histograms have been 
constructed using 60 random realizations each. The sohd black line is for a model with no 
scatter, the red dotted line for a model with a uniform dispersion in ki and the green dashed 
line for a model with a lognormal distribution (with a^: = 1). 



^nt^' distribution from Monte Carlo realizations with a dust-to-gas ratio scatter 
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Fig. 7. — Like fig. 6, but for syntlietic observations witfi a lognormal ki distribution (a* = 1) 
and a large correlation (7 = ±0.8) between ki and Nhi (green dashed and red dotted lines) 
vs. realizations with no scatter (black solid line). 
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Fig. 8. — Cumulative neutral gas density in DLA systems versus column density for the SDSS 
DLA DR_3 {observed dotted, intrinsic with = 1 solid). Superimposed to the intrinsic 
curve we show the profile from the best fitting gamma function /j (short dashed) and from 
a power law fi (long dashed) with slope —1.95: the discrete evaluation for the intrinsic 
distribution does not extend to sufficiently high column densities to allow convergence. 
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Table 1: Maximum likelihood parameters for gamma function fi 



OLk 


^1 




mi 




log{C) 






0.000 


3.03 10- 


-3 


-1.80 


3.00 


0.00 


0.821^5 10- 


-3 


0.125 


2.61 10- 


-3 


-1.80 


3.23 


0.16 


0.851^5 10- 


-3 


0.250 


2.49 10" 


-3 


-1.78 


3.40 


0.32 


0.891°:°^ 10- 


-3 


0.375 


2.21 10" 


-3 


-1.78 


3.61 


0.46 


10- 


-3 


0.500 


1.70 10" 


-3 


-1.79 


4.15 


0.66 


0.97+°:°« 10- 


-3 


0.625 


1.47 10" 


-3 


-1.79 


4.51 


0.83 


i.oi;°:J° 10- 


-3 


0.750 


1.10 10" 


-3 


-1.79 


5.27 


1.00 


1.08+°:}^ 10- 


-3 


0.875 


1.02 10" 


-3 


-1.78 


5.60 


1.18 


1.13l}3 10- 


-3 


1.000 


6.24 10- 


-4 


-1.80 


7.22 


1.33 


l.2ht',^ 10- 


-3 



Note. — Summary table with the maximum hkehhood parameters for the intrinsic distribution of DLA 
systems fitted with a gamma function and using different values for a^. The parameters rji and are given 
adopting units of 10^^ cm~^ for Nhi- The likelihood C is normalized to the value of the best fitting model 
(gamma function) with no dust. 
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Table 2: Maximum likelihood parameters for a power law function /j 







V2 




m2 


log{£) 


0.000 


1 


27 


10- 


-2 


-2.19 


-9.90 


0.125 


1 


36 


10- 


-2 


-2.14 


-7.07 


0.250 


1 


43 


10- 


-2 


-2.10 


-5.26 


0.375 


1 


49 


10- 


-2 


-2.07 


-3.84 


0.500 


1 


54 


10- 


-2 


-2.05 


-2.70 


0.625 


1 


61 


10- 


-2 


-2.01 


-1.73 


0.750 


1 


67 


10- 


-2 


-1.99 


-1.95 


0.875 


1 


74 


10- 


-2 


-1.96 


-0.24 


1.000 


1 


77 


10- 


-2 


-1.95 


0.29 


1.250 


1 


91 


10- 


-2 


-1.90 


1.19 


1.500 


2 


05 


10- 


-2 


-1.85 


1.82 


2.000 


2 


35 


10- 


-2 


-1.75 


2.39 


4.000 


3 


75 


10- 


-2 


-1.45 


-1.12 



Note. — Summary table with the maximum hkehhood parameters for the intrinsic distribution of DLA 
systems fitted with a power law function and using different values for a^. The parameters rji and are 
given adopting units of lO^^cm"^ for Nhi- The likelihood C is normalized to the value of the best fitting 
model (gamma function) with no dust. 
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Table 3: Maximum likelihood parameters for the CORALS survey 



6.33 10- 


-6 


mi 
-2.13 


39.09 


log{C) 



^(DLA) 

1.2 10-=^ 


1.87 10- 


-2 


1712 

-1.80 




log{£) 
-1.15 




1.47 10- 


-2 


-2.20 




-0.05 




1.11 10- 


-2 


-2.50 




-0.55 





Note. — Summary table with the maximum hkehhood parameters for the CORALS survey modehng. 
The Ukelihood value is normalized to the value given by the best fitting CORALS model. Other units are 
hke in Table 1 
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Table 4: Synthetic observations with scatter in the dust-to-gas ratio 



Model 






10^^(1]^/^^) 


io-V(Qa/) 


IO'^AQh/ 


Uniform 50%-150% 




LOO 


0.11 


0.03 


LogNor (7* = 0. 


5 




0.97 


0.08 


0.00 


LogNor cr* = 1 






0.99 


0.09 


0.02 


LogNor (7,^ = 1. 


5 




L03 


0.10 


0.05 


LogNor (7* = 1; 


7 = 


+0.4 


0.97 


0.10 


0.00 


LogNor (7* = 1; 


7 = 


-0.4 


1.05 


0.10 


0.08 


LogNor cr^, = 1; 


7 = 


+0.8 


0.94 


0.07 


-0.03 


LogNor (7* = 1; 


7 = 


-0.8 


1.11 


0.13 


0.14 



Note. — Summary of the synthetic observations with a scatter in the dust-to-gas ratio ki. The underlying 
column density distribution is the best fitting SDSS gamma function when = 0.5 (A^-y = 4.15 • lO^^cm^^, 
mi = —1.79, il^^'^-' = 0.97 • lO^''). The first series of models employs a uncorrelated dispersion in ki 
(uniform and lognormal, see Sec. 6), the second series also assumes a correlation between ki and Nhj. For 
each model we have analyzed 60 random realizations to recover O^^"*' reporting the mean value obtained 
from the maximum likelihood analysis with average ki{z) (second column), the standard deviation (third 
column) and the bias, i.e. the average value minus the nominal model value (fourth column). 
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Table 5: Discrete evaluation of fl 







^max 


dXo 


Ndla 










2.2 


5.5 


7333.1 


525 


1.04 ±0.08 


0.91 ±0.07 


0.82 ±0.06 




2.2 


2.5 


1589.2 


83 


0.53 ±0.07 


0.48 ±0.06 


0.45 ± 0.06 


SDSS: 


2.5 


3.0 


2697.6 


189 


1.04 ±0.14 


0.90 ±0.11 


0.79 ± 0.09 




3.0 


3.5 


1836.3 


152 


1.27 ±0.18 


1.13±0.15 


1.01 ±0.13 




3.5 


4.0 


877.3 


69 


1.47 ±0.30 


1.27 ±0.24 


1.11 ±0.19 




4.0 


5.5 


332.6 


33 


1.13 ±0.31 


1.06 ±0.27 


1.00 ±0.24 






^max 


dXo 


Ndla 


lo^nl^,^^^ 






CORALS: 


1.8 


3.5 


200.8 


17 


1.20 ±0.50 







Note. — ^Hi evaluated in the discrete limit for different rcdshift intervals from the SDSS and CORALS 
data. The first two columns give the rcdshift interval considered, the third the corresponding pathlcngth 
probed by the data {dXo)- Ndla is the number of systems identified in the interval, fi^^'^'' the intrinsic 
density of neutral gas (evaluated using = 1, = 0.5 and = 0). The uncertainty is evaluated with a 
'Jack knife' analysis. 



